#ifndef QUANTIZE_H
#define QUANTIZE_H
#include <fcntl.h>
#include <iostream>
#include <sys/mman.h>
#include <unistd.h>
#include <boost/numeric/mtl/mtl.hpp>
#include <cfloat>
#include <flann.h>
#include "configure.h"
#include "myutility.h"

const float EPS = 1e-6;
using namespace mtl::matrix;
class Quantizer
{
    private:
    float* clust;
    dense2D<float> hamatrix;
    FLANN_INDEX findex;
    FLANNParameters fp;
    int topnn;

    float* hashclust;

    public:
    // buffer for format a single siftgeo file
    int* nn;
    float* dist;
    float* hashsift;
    char* hamcode;

    Quantizer()
    {

    }
    Quantizer(string clust_path,string hamatrix_path,int _topnn = 1,int maxrecords = MAXRECORDS_PER_SIFTGEO)
    {
        cout << "quantizer constructor " << endl;

        clust = NULL;
        load_mat(clust_path,clust);
        load_hamatrix(hamatrix_path,hamatrix);
        hashclust = NULL;


        topnn = _topnn;
        nn = new int[maxrecords*topnn];
        dist = new float[maxrecords*topnn];
        hashsift = new float[maxrecords*HAMMING_DIM];
        hamcode = new char[maxrecords*(HAMMING_DIM>>3)];

        float speedup;
        fp.algorithm = LINEAR;
        fp.target_precision = 1.0;
        fp.trees = 8;
        fp.memory_weight = 0.01;
        fp.build_weight = 0.01;
        fp.log_destination = NULL;
        fp.log_level = LOG_ERROR;
        findex = flann_build_index(clust,LEXICON_SIZE,SIFT_DIM,&speedup,&fp);
    }
    Quantizer(string clust_path,string hamatrix_path,string hashclust_path,int _topnn = 1,int maxrecords = MAXRECORDS_PER_SIFTGEO)
    {
        cout << "quantizer constructor " << endl;

        clust = NULL;
        load_mat(clust_path,clust);
        load_hamatrix(hamatrix_path,hamatrix);
        hashclust = NULL;


        topnn = _topnn;
        nn = new int[maxrecords*topnn];
        dist = new float[maxrecords*topnn];
        hashsift = new float[maxrecords*HAMMING_DIM];
        hamcode = new char[maxrecords*(HAMMING_DIM>>3)];

        float speedup;
        fp.algorithm = LINEAR;
        fp.target_precision = 1.0;
        fp.trees = 8;
        fp.memory_weight = 0.01;
        fp.build_weight = 0.01;
        fp.log_destination = NULL;
        fp.log_level = LOG_ERROR;
        findex = flann_build_index(clust,LEXICON_SIZE,SIFT_DIM,&speedup,&fp);
        load_mat(hashclust_path,hashclust);
    }
    ~Quantizer()
    {
        cout << "quantizer destructor " << endl;
    }
    void free()
    {
        if(clust != NULL) delete clust;
        if(hashclust != NULL) delete hashclust;
        if(nn != NULL)  delete[] nn;
        if(dist!= NULL )delete[] dist;
        if(hashsift != NULL)delete[] hashsift;
        if(hamcode != NULL)delete[] hamcode;
        flann_free_index(findex,&fp);
    }
    static void load_mat(string matpath,float* &mat)
    {
        int row,col;
        struct stat statbuf;
        stat(matpath.c_str(),&statbuf);
        int bytenum = statbuf.st_size;
        int fid = open(matpath.c_str(),O_RDONLY);
        char *back_buf = reinterpret_cast<char*>(mmap(NULL,bytenum,PROT_READ,MAP_PRIVATE,fid,0));
        row = *(reinterpret_cast<int*>(back_buf));
        col = *(reinterpret_cast<int*>(back_buf)+1);
        cout << "row: " << row << " col: " << col << endl;
        float* buf = reinterpret_cast<float*>(back_buf)+2;
        mat = new float[row*col];
        memcpy(mat,buf,sizeof(float)*row*col);
        munmap(back_buf,bytenum);
        close(fid);
    }
    static void load_hamatrix(string hamatrix_path,dense2D<float> &hamatrix)
    {
        int row,col;
        int fid = open(hamatrix_path.c_str(),O_RDONLY);
        struct stat statbuf;
        stat(hamatrix_path.c_str(),&statbuf);
        float *buf = reinterpret_cast<float*>(mmap(NULL,statbuf.st_size,PROT_READ,MAP_PRIVATE,fid,0));
        char* back_buf = reinterpret_cast<char*>(buf);
        row = *(reinterpret_cast<int*>(buf));
        ++buf;
        col = *(reinterpret_cast<int*>(buf));
        ++buf;
        cout << "row: " << row << " col: " << col << endl;
        hamatrix = dense2D<float>(row,col);
        for(int i=0; i<row; ++i)
            for(int j=0; j<col; ++j)
                hamatrix[i][j] = *(buf++);
        munmap(back_buf,statbuf.st_size);
        close(fid);
       // is_orthogonal(hamatrix);
    }
    void findnn(float* vec,int vecnum)
    {
        cout << " use flann library :) " << endl;
        flann_find_nearest_neighbors_index(findex,vec,vecnum,nn,dist,topnn,fp.checks,&fp);
    }
    void _findnn(float* vec,int vecnum)
    {
        cout << " use my poor implemention :( " << endl;
        float dis = 0.0;
        float mindis = FLT_MAX;
        float *vp, *cp,*pp;
        vp = vec;
        print_matrix(clust,LEXICON_SIZE,SIFT_DIM);
        for(int i=0; i<vecnum; ++i,vp+=SIFT_DIM)
        {
            mindis = FLT_MAX;
            cp = clust;
            for(int j=0; j<LEXICON_SIZE; ++j)
            {
                dis = 0;
                pp = vp;
                for(int k=0; k<SIFT_DIM; ++k,++cp,++pp)
                dis += (*pp-*cp)*(*pp-*cp);
                if(dis < mindis)
                {
                    mindis = dis;
                    nn[i] = j;
                }
            }
        }
    }
    void gen_hashsift(float* sift,int siftnum)
    {
       dense2D<float> siftmatrix(siftnum,SIFT_DIM);
        for(int i=0,k=0; i<siftnum; ++i)
            for(int j=0; j<SIFT_DIM; ++j,++k)
                siftmatrix[i][j] = sift[k];
        dense2D<float> hashsiftmatrix(siftnum,HAMMING_DIM);
        hashsiftmatrix = siftmatrix * hamatrix;
        for(int i=0,k=0; i<siftnum; ++i)
            for(int j=0; j<HAMMING_DIM; ++j,++k)
                hashsift[k] = hashsiftmatrix[i][j];
    }
    void gen_hamcode(float* hashsift, int num)
    {
        int bytes = HAMMING_DIM >>3;
        char mark = 0;
        float* hclustptr = NULL;
        float* hsiftptr = hashsift;
        char* hamcodeptr = hamcode;
        for(int i=0; i<num; ++i)
        {
            hclustptr = hashclust + HAMMING_DIM * nn[i];
            for(int j=0; j<bytes; ++j)
            {
                mark = 0;
                for(int k=0; k<8; ++k,++hclustptr,++hsiftptr)
                {
                    if(*hsiftptr >= *hclustptr)
                        mark |= (1<<k);
                }
                *(hamcodeptr++) = mark;
            }
        }
    }
    bool is_orthogonal(dense2D<float> &hamatrix)
    {
        int row = num_rows(hamatrix);
        int col = num_cols(hamatrix);
        dense2D<float> thamatrix(col,row);
        for(int i=0; i<row; ++i)
            for(int j=0; j<col; ++j)
                thamatrix[j][i] = hamatrix[i][j];
        dense2D<float> mul(col,col);
        mul = thamatrix * hamatrix;
        for(int i=0; i<col; ++i)
            {
                for(int j=0; j<col; ++j)
                {
                    if(i==j)
                    {
                        if(fabs(mul[i][j] -1) > EPS)return false;
                    }
                    else if(fabs(mul[i][j])>EPS) return false;
                }
            }
        cout << "is a orthogonal matrix" << endl;
        return true;
    }
};
#endif
